Ab initio description of nonlinear dynamics of coupled microdisk resonators with 

application to self-trapping dynamics 



O 

U 

Oh' 
< 
(N 



O 

H— > . 

Or 

q 
o 

43 ■ 



> 

o 

(N 
(N 

O 



X 



Hamidreza Ramezani 1 , Tsampikos Kottos 1,3 , V. Shuvayev 2 , and L. Deych 2 
1 Department of Physics, Wesleyan University, Middletown, Connecticut 06459, USA 
2 Department of Physics, Queens College of the City University of New York (CUNY) Flushing, NY 11367 and 
Max-Planck-Institut fur Dynamik und Selbstorganisation, Gottingen, Germany 
(Dated: January 19, 2013) 

Ab initio approach is used to describe the time evolution of the amplitudes of whispering gallery 
modes in a system of coupled microdisk resonators with Kerr nonlinearity. It is shown that this 
system demonstrates a transition between Josephson-like nonlinear oscillations and self-trapping 
behavior. Manifestation of this transition in the dynamics of radiative losses is studied. 



I. INTRODUCTION 

Whispering gallery modes (WGM) are optical excita- 
tions that exist in axially symmetrical optical resonators. 
In the geometric optic framework they can be described 
as waves propagating at large incidence angles along the 
surface of the resonator and trapped in an optical po- 
tential arising due to total internal reflection. Due to 
high Q-factors and small mode volumes of WGMs, one 
should expect a significant enhancement of nonlinear ef- 
fects in such structures [lH3|. Indeed, a variety of nonlin- 
ear phenomena such as ultra-low threshold Raman las- 
ing [|| , parametrical optical oscillations Q, and optical 
comb generation [y] have been experimentally demon- 
strated in WGM resonators. Theoretical understanding 
of these effects, however, relied mostly either on purely 
phenomenological approaches similar to that of Ref . 7J > 
or employed approximations of the mean-field type [!, 8| . 

When several WGM resonators are placed in the 
proximity of one another, they become optically cou- 
pled due to overlapping evanescent tails of their respec- 
tive modes. The resulting collective optical excitations 
in linear regime have been observed experimentally in 
coupled microspheres and micro-disks 0-[Hl and de- 
scribed theoretically within the framework of ab initio 
approaches pH4l^ |. At the same time, the nonlinear 
optics of coupled microresonators is still in its infancy. 
While lately there was a significant uptick in theoreti- 
cal attention to this field |15l4l8l |. no experiments in this 
area have yet been reported. This rise of interest to op- 
tically coupled resonators is fueled by expectations that 
these systems would allow for a convenient experimen- 
tal studies by optical means of some important nonlinear 
effects of generic nature. In this work we demonstrate 
theoretically that coupled microdisk resonators with in- 
stantaneous Kerr nonlinearity is a convenient system, in 
which such nonlinear effects as Josephson oscillations and 
self-trapping can be demonstrated. This system is signif- 
icantly simpler for ex per imental realization than those 
considered in Ref. [TBI. Il7l. Il8j. and we expect, therefore, 
that this work will motivate experimental studies of non- 
linear optical properties of microdisks and other coupled 
WGM resonators. 

Unlike previous theoretical works relying on phe- 



nomenological models |15l - ll8l |. we use an ab initio ap- 
proach to derive equations describing nonlinear dynam- 
ics of collective WGM excitations in this system. We 
show that under certain conditions these equations can 
be presented in the form of modified discrete nonlinear 
Schrodinger equation (DNLSE). DNLSE is a prominent 
model known to describe successfully the dynamics of 
systems as diverse as Bose-Einstein Condensates (BEC) 
in optical lattices, polarons, optical fibers, and biologi- 
cal molecules [17H26l | . DNLSE describes such nonlinear 
phenomena as self-trapping, discrete breathers, etc. and 
demonstrate anomalously slow relaxation dynamics |27| 
(due to boundary losses) similar to the one appeared in 
glass materials. This relaxation dynamics has generated 
recently a great deal of theoretical interest in the context 
of cold atoms [28l - l3l| and coupled cavities fl7l [l8j . At 
the same time, while self-trapping in systems with con- 
serving number of particles or energy has been observed 
experimentally 24] , it is still yet to be observed in in any 
physical realization of the DNLSE with losses. 

In this work, we demonstrate that the system of cou- 
pled microdisk cavities provides a convenient platform 
for studying nonlinear dynamics in the presence of re- 
laxation. While radiative losses in such a structure are 
intrinsic to each disk, one can introduce an imbalance of 
the losses required for anomalous relaxation by coupling 
one of the disks to a tapered fiber. Such an arrangement, 
which is usually used to excite WGMs, is known to reduce 
the Q-factor of the fiber-coupled resonator 32] providing 
thereby a faster relaxation channel through the boundary 
of the structure. We demonstrate that at certain thresh- 
old values of initial conditions/nonlinearity the dynam- 
ics of the system exhibits a transition between nonlinear 
Josephson-like oscillations of light intensity between the 
disks and a self-trapping behavior. We also show that 
the transition between these two regimes is manifested 
in different kinetics of radiative losses in Josephson and 
self-trapping regimes. 



II. AB INITIO DYNAMICAL EQUATIONS 

While we focus in this paper on a nonlinear dynamic 
of two laterally coupled microdisk resonators, the general 
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dynamic equations can be derived for a more general case 
of N disks, provided they all lie in the same plane, which 
is parallel to their surfaces. In the spectral region of 
high-Q WGM resonances electromagnetic field of a sin- 
gle disk can be described by a simplified two-dimensional 
model [33] , in which its electric (for TM modes) or mag- 
netic (for TE modes) field is perpendicular to disk's sur- 
face and is, therefore, characterized by a single function 
F(*)(r, t). In this approximation, disks' refractive index is 
replaced by an effective parameter, rid, determined from a 
self-consistency condition [33| . Two-dimensional approx- 
imation was also used to describe electromagnetic field 
of the planar system of several disks (optical molecules) 
using coupled integral equations [13, [S3 as well as modal 
expansion approach 11]. While the formulation in terms 
of integral equations is convenient for numerical simu- 
lations of planar structures, especially with non-circular 
elements, the modal expansion is more appropriate for 
developing ah initio theoretical description of nonlinear 
effects in WGM resonators. 



defined by equations 
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L A^G(r - r') = -5(r - r'), r < R 
f (2) 



V 2 G(r - r') + — G(r - r') = -S(r - r'), r>R 

and boundary conditions 

G(R-0,r') = G(R + 0,r') 



dG(r, r 1 ) 



dG(r, r') 



® r r=R~0 

d In [r^ 2 G(r, r')] 



r=R+0 



(3) 



di 



where the last of these expressions establishes outgoing 
field condition at infinity. Using Green's theorem [36| one 
can express the field in the interior of the disk in terms 
of the Green's function defined in Eq. [2] and boundary 
values of the field and its derivative: 



Single disk in presence of external incident field 
and internal polarization 



Before developing the theory of the nonlinear dynamics 
of WGMs in multiple disks, one has to consider the case 
of a single disk. We will assume that the nonlinearity 
does not couple TM and TE polarizations and limit our 
consideration to the TM case. This assumption is valid 
for disks made of materials with small optical anisotropy 
and in the absence wave mixing processes. In the spectral 
domain two-dimensional wave equation for the Fourier 
transformed electric field can be written as 



V 2 F(r,0, W ) + ^F(r,0, W ) = 



AlTUJ 1 



-P,r <R 



UJ 



(1) 



V 2 F(r, cj>, oj) + -irF(r, 4>, u) = 0,r>R 



where r, and <p are, respectively, radial and angular po- 
lar coordinates, defined in a coordinate system with the 
origin at the disk's center, to is a spectral parameter, not 
yet identified with any physical frequency, is the ef- 
fective refractive index of the disk, c is vacuum speed of 
light, and P is the nonlinear polarization. Usually, Eq[T] 
is solved in the context of the scattering problem, when 
the boundary conditions are of inhomogeneous type de- 
fined in the presence of an external "incident" field. At 
the same time, nonlinear effects are usually described in 
terms of dynamics of " modal amplitudes" , with modes 
being solutions of a linear problem with homogeneous 
(no incident field) boundary conditions. The difficulty 
for the nonlinear theory of WGMs, which are usually ex- 
cited by an "incident" field, is the necessity to reconcile 
these two different types of problem. In order to achieve 
this we, following Ref. [36|, introduce a Green's function 



F m {r) = 
- 2-kR 
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G m (r',r)P m (r r ydr> 



F m (Rf-^-G M dF ^ 



dr 1 



dr' 



r'=R 

(4) 



where all quantities with subindex m are angular Fourier 
coefficients of their respective functions: < • > m = 
(l/27r) < ■ > ((f)) exp [im(j}]d4>. The values of the 
field and its derivative at the boundary of the disk, 
r = R, are defined via Maxwell boundary conditions 
demanding continuity of F and its derivative with re- 
spect to the radial variable. In the presence of the in- 
cident field, F inc (r,(f)), the field outside of the disk is 
Finc(r, (f>) + F sc (r, <£>), where F sc (r, (f>) is the field scattered 
by the disk. Taking into account standard conditions for 
the incident field to be finite at r = and for the scat- 
tered field to be outgoing at infinity, these fields can be 
presented as linear combinations of the first order Bessel 
(J m ) and Hankel (H m ) functions: 

Fine = ^ &'mJm(kr S ) exp (im(j)) (5) 

7) 1 

F sc = y^b m H m (kr)exTp(im(f>) (6) 



where fc = lo/c. Combining Eq. U with the boundary 
conditions as well as with Eq. [5] and [S] one arrives at the 
following expression for the internal (r < R) field of the 
disk: 



Fm(p) 



8tt 2 x 2 



2i 



G m {p\p)P m ( P ')p'dp' 
a m J m (ridxp) 



(7) 
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where we introduced dimensionless spectral parameter 
x = kR and radial coordinate p — r/R. Function £ m (x) 
in the second term of Eq. [7] is defined as 

£ m (x) = n d H m (x)J' m (n d x) - J m (ji d x)H' m (x) (8) 

where f'(z) = df/dz. 

In the absence of the polarization source, Eq. [7] re- 
produces a standard solution for the internal field in the 
linear single disk scattering problem. The scattered field 
in this case is also found in a standard form 



where a m is a linear scattering amplitude defined as 



(9) 



£m(x) 



(10) 



where p m (x) is 

Pm(x) = J' m {x)J m {n d x) - n d j' m (n d x)J m (x). (11) 
The poles of a m (x), which we will designate as Xm]s — 

(r) 

i^fm's, are given by complex- valued solutions of equation 
£(x) = and define spectral positions, Xm,s, and widths, 

(r) 

jm's, of linear WGM resonances in a single disk. The sub- 
index s here distinguishes resonances with different radial 
distributions of the respective internal fields: the value of 
this index determines the number of maxima of the field 
in the radial direction. The scattering amplitude has an 

important property: a m (xm]s) = — 1, which results in the 
following single pole approximation valid in the vicinity 
of a selected resonance: 



■ 0) 



Jr) 



, ■ 0) 



(12) 



Eq.[9land[T0lsolve the linear scattering problem and de- 
scribe resonant response of the disk to the external exci- 

(r) 

tation. However, complex frequencies x m \s are not eigen- 
frequencies of the disk, and their respective field distri- 
butions are not its eigenfunctions. In order to define true 
normal modes of the system together with their eigenval- 
ues, one has to solve the wave equation with boundary 
conditions formulated in the absence of the incident field 
(homogeneous boundary conditions). The open nature 
of the optical resonators makes this problem nontrivial. 
Several different approaches have been developed to in- 
troduce a system of functions, generalizing the concept 
of modes, that could be used as a basis to describe open 
systems (see, for instance, recent review [37j). For our 
goals, the most convenient is the approach based on so 
called Constant Flux Modes, re-introduced in the opti- 
cal context in Ref. 38] and used recently in Ref. [39| . 
These functions are defined as solutions of the following 
equations 

V 2 4, s (p, 0) + n 2 d x 2 m ^ m ^(r, (p)=0r <R 

72„/. /„ A\ , „2v ' ^ " - " ' J 



with the same boundary conditions as those introduced 
in Eq. [3] for the Green's function. It is crucial that out- 
side of the disk these functions obey the wave equation 
with generic spectral parameter x rather than with the 
eigenvalue x mjS . Only functions defined this way form, 
in combination with their adjoint counterparts, a com- 
plete bi-orthogonal set with an inner product and norm 
defined over the interior of the disk: 



ipm, s (p) \tfi, p {p)] * Pdp = N? n s 5 m .i5p, s (14) 



Adjoint functions ifi^ p (p) appearing in Eq. [14] are solu- 
tions of equations 

V 2 ^ m Jp,<P)+nl[x 2 m J*^ s (r, ( l ) )=Or<R 
V 2 i, s (^)+i 2 i,(^)=0,r>iJ 

with incoming boundary conditions at infinity. The 
eigenvalues x m . s are defined by equation 

^dXm^s H m (x) J m {jidXjn^s) xJ m (n d x rns )H m (x) 

and are complex- valued functions of the spectral parame- 
ter x, which becomes an integration variable upon trans- 
formation of the fields back to the time-domain. 

Since eigenfunctions ip miS obey the same boundary 
conditions as the Green's function, G m (p,p r ), they can 
be used to generate its spectral expansion 

„ , R 2 ^rn,s(x m . s n. d p)tp* m s (Xm, S n d p') 

G m {p,p) = -t^T2 2^ 



2TTU 2 



2 2 

x — x m 



The same set of functions can be used to present the 
internal field of the disk 

F m (p, x) = D mtS (x)ip mtS (x m<s ndp) (18) 



yielding, in combination with Eq. 01 Eq. [T7] and bi- 
orthogonality relation, Eq. [14l the following system of 
equations for modal amplitudes D m s : 



D m ,s (x x m _) 



4?rx 2 An,s(x) 
ni 



p m {x) 

Here we introduced projections of the nonlinear polariza- 
tion, P m ,s, and incident field, A m ^ s , onto eigenfunctions 
of the disk: 



Pjn,s (x) 



P m (x,p)tp (x miS n d p)pdp (20) 



V ip m>s (r,4>) +x ip m , s (r, 



Q,r>R 



A m ,s(x) = — / J m (n d xp)ip m s (x m , s n d p) pdp (21) 



(22) 



and a parameter K m ,s defined as 
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Parameter K m>s in the last term of Eq. [19] characterizes 
the efficiency of coupling of the incident radiation to the 
internal modes of the disk. The structure of this param- 
eter, which we will call external coupling parameter, re- 
flects distinction between external and internal excitation 
mechanisms: external incident field induces response at 
scattering resonances, while internal excitation (polar- 
ization term) generates response at the eigenfrequencies. 
Taking into account Eq.[12]one can present K m . s in vicin- 
ity of x m , s as 



2^7, 



m.s^rn.s 



(x- 



.) 



Jr) 



■ (r) 



(23) 



For resonances with low Q-factors the difference between 
eigenfrequencies and scattering resonances can be quite 
significant resulting in a strong frequency dependence of 
the coupling to the external field. However, for high- 
Q resonances, numerical calculations indicate (see also 

Ref. [11|) that lZex mtS — Xm]s as well as 2mx„ ltS — 7™! 
are much smaller than the individual imaginary parts of 
these quantities. In this case k can be simplified to the 
following form 



-2ij 



(24) 



which shows that efficiency of coupling to the modes of 
the disk from outside is proportional to the width of 
these modes. This physically clear result expresses the 
fact that one cannot couple from outside to a completely 
closed system. However, to the best of our knowledge, 
this result has not been derived previously from "first 
principles" and in many cases (see for instance, Ref. (ipj ) 
a factor representing the width of the resonances was in- 
troduced into the respective term "by hands". 



B. Nonlinear dynamics of multiple optically 
coupled disks in the frequency domain 

1. General equations 

The case of multiple disks can be treated by applying 
Eq. 033 to each disk of the structure and including in the 



term containing the incident field contributions from the 
fields scattered by all other disks [111 ]: 



v^p n 



(25) 



Here indexes p,v enumerate disks comprising the struc- 
ture, and t v m p _ n = e^ n - m ^-"H m - n (xR V:P /R) describes 
optical coupling between the disks, where R VjP and 9 V p 
are, respectively, radial and angular polar coordinates of 
u-th disk relative to the p-th one. This terms arises from 
transformation of the field scattered by v-th disk to the 
coordinate system centered at p-th disk using Graf's ad- 
dition theorem for Hankel functions [4l[ . If all the disks 
are arranged along a single line, which is chosen as a po- 
lar axis of the coordinate system, one has 6 pv = ir for 
p < v and 9 pv — for p > v. 

Eq.[l9]in the multiple disk case must be complimented 
by an equation relating scattering coefficients \P n of dif- 
ferent disks to each other. In the absence of nonlinearity 
this equation is again obtained by substituting Eq. 1251 to 
the single disk Eq. |9] which gives [llj 



K, = ot m {x) 



n v^p 



(26) 



While the presence of nonlinear polarization in Eq. Q] 
modifies Eq. 1261 this modification can be neglected since 
it is proportional to two relatively weak effects: inter- 
disk coupling and single-disk nonlinearity Only when 
frequency shifts due to any of these effects becomes com- 
parable with linear frequencies of the resonator, nonlinear 
corrections to Eq.[2S] should be taken into account. Com- 
bining Eq. Q1J] with Eq. [23] and Eq. [211 one can eliminate 
coefficients of the scattered field, 6™ from the equation 
for the modal amplitudes and obtain the closed equation 
for the latter: 



Dm,s K m,s ( x ) y 1 y 1 ^m—n^n, 



(X) 



Aitx 2 p r , 



_pp 



(27) 



where we introduced renormalized modal amplitude according to 
-D^j s related to the original amplitude defined in Eq. [18] 

I) 1 ' 

A. 



= -P-Btn,. (28) 



This amplitude describes the field outside of the res- 
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onator produced by an internal mode with amplitude 
DJ^ s . It is important to note that the external coupling 
coefficient n m ,s appears in Eq. [27] not only in the term 
containing incident field, but also in the one describing 
inter-disk coupling. This is a reflection of the obvious, 
but still not always appreciated, fact that the optical 
coupling between disks occurs via scattered rather than 
internal fields. 

Equation [27] must be, of course, complimented by 
an expression for the nonlinear polarization in terms of 
modal amplitudes, D m s . Assuming nonlinearity of Kerr 



type, we can present nonlinear polarization term P(ui) as 

cn 2 d n 2 fdx 1 dx 2 , „ . , . , , 
16tt z J 2n 2n 

where x% t 2 = x—x\ +x 2 and we used approach of Ref. [42| 
to express Kerr nonlinear polarization in terms of stan- 
dard nonlinear refractive index n 2 . Combining Eg 1291 
with modal expansion Eq. [18] one can finally derive the 
following expression for the component of the nonlinear 
polarization P lritS in terms of modal amplitudes D m , a (x): 



5 (» 



m-m2+m 3 ,mi J Vi,i/2,i'3 



dx\ dx 2 
~27T~27r 



A/i (^1,2) A 2 (xi)Dl (x 2 ) 



(30) 



where for shortness we combined double indexes mj,s, coupling coefficients A™'* 9 „ 3 defined as 
in single indexes Vi, and introduced nonlinear inter-mode 



cn d n 2 



A %',t 2 ,»/3 = 16?r 2 / d PPi J i'i( x ^i n dp)i'u 2 (x^ 2 n d p)i; U3 (x U3 n d p)Tp ms (x m ^n d p) (31) 



The derived equations [27] - [31] provide an ab initio foun- 
dation for studying a variety of nonlinear processes in the 
the system of coupled disk resonators with Kerr nonlin- 
earity. The next step in development of the theory would 
include Fourier transformation of the derived equations 
back to time-domain. This procedure cannot be car- 
ried out exactly and requires approximations of a "slow 
changing amplitude" type. The implementation of the 
latter, however, depends on the type of nonlinear pro- 
cesses being studied. In this paper we are interested in 
self-trapping dynamics in the system of coupled disks, 
thus in the subsequent sections of the paper we will ap- 
ply Eq. [57] - [3T] to this particular problem. 



ynamic of modal amplitudes in a double-disk optical 
"molecule" 



In what follows we limit our consideration to the case 
of only two disks, which we will study in the so called 
resonant approximation [Tl| . This approximation im- 
plies that we only take into account coupling between 
degenerate phase-matched modes of the disks. These are 
clockwise and counterclockwise modes characterized by 
azimuthal numbers of opposite signs. The phase match- 
ing means that mode with m = M in disk 1 couples 



only to the mode m = -M of the disk 2 and vice versa. 
On formal level this is justified by inequality t 2 M 3> to 
valid for M 3> 1. In this approximation Eq . |2"7I takes the 
following form: 



(a; 2 — x 2 M ) D\ 4 — K„ liS (x)t2 M D'L M (x) = 
[x 2 — x 2 M ~} D 2 M — K mtS (x)t 2 ' M D^ M (x) = 
(x 2 — x\j) D 1 _ M — K m ^ s (x)t2 M b 2 M {x) = 
(x 2 — x 2 M ) D 2 i M — K m , s (x)t 2 ' M b\ t {x) = 



4ttx p M 51 
n d A M 
4ttx 2 p m ~ 2 



M 

n d A M 
4ttx 2 Pm pi 

n 2 d A M 
4irx 2 p M „ 



—M 



n 2 A M 



P; 



-M 



(32) 



where we abridged notations by replacing double index 
m, s with a single index M having in mind that in the res- 
onant approximation all involved modes have the same 
radial index set to be equal to unity. In Eq. [35] we also 
omitted the term with the incident field since we model 
excitations of the dynamics in this paper by initial con- 
ditions rather than by an external field. 

Transformation of Eq. [35] to the time domain is based 
on a slow amplitude approximation, when the nonlinear 
dynamics is presented as a slow modulation of fast linear 
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oscillations of the field: 



1 



D p M (t) = ^(S P M (t)e- tXot +c.c.) 



(33) 



where S P j is a slow changing amplitude, xq is the fre- 
quency of the linear dynamics, and c.c. means "complex 
conjugated" term. In the spectral domain this approxi- 
mation has the following form 



D p M (x) = - [S p M {x)S(x - x Q ) + [S p M (x)}* 8(x + x )] 



The choice of Xq depends on the problem at hands. For 
instance, in the presence of an external incident field, 
xq would be set by its frequency, while in the case of 
initial value problem, considered in this work, the lin- 
ear dynamics is determined by the poles of the Green's 
function, Eq. [T7] These poles are solutions of equation 
x m ,s(x) = x, and it is clear from Eq. [8] and Eq. [TBI that 
they coincide with the scattering resonances. Thus, in 
this case x = lZe[x$] and Eq. [M]for the external cou- 
pling parameter becomes exact. 

When using Eq. [M] to convert system [32] to the time 
domain, one needs to take into account the dependence 
of the eigenfrequencies on the spectral parameter x and 
expand it in the power series around x = xq ■ In the slow 
amplitude approximation one only keeps term linear in 
x, which correspond to the first order time derivative in 
the time-domain: 



X M (x) ~ Xq + (x - Xq) 



dx 



A I 



dx 



While the dxiw/dx term affects the coefficient in front 
of term with the first order time derivative of the modal 
amplitude, its numerical estimates showed that in the 
case under consideration in this work it can be neglected. 

Time-domain expression for the nonlinear polarization 
terms is obtained by Fourier-transforming Eq. [30] with 
the help of Eq . [3~T1 and restricting summation over modes 
only to those with s = 1 and m = ±M. In doing so, 
we also neglect all fast oscillating terms at combination 
frequencies. Combining the result with the Fourier trans- 
form of the linear part of Eq[321we obtain the final set of 
time-dependent equations for modal amplitudes 



■ ^±m( T ) . ■ a (v) f \ . , a(p) , \ 

i j- 1- «7m^±m(t) + h M S^{r) 



-C,T M S± M 



l4tl 2 



215. 



(v) 1 2 

=fmI 



(35) 



In these equations we introduced dimensionless time r = 
tc/R and linear coupling coefficient hu = i^lQ t 2 Af (recall 

(12) 

that the imaginary part of Hankel function entering t 2M 
is much larger than its real part for M > 1, so that 
hyi has very small imaginary part). Kerr nonlinearity 
results in self-coupling of modes and cross - coupling be- 
tween counter-propagating modes. The strength of non- 
linear interactions is characterized by material parame- 
ter £ = (cn2Xo)/(47r), and the dimensionless resonator 
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FIG. 1: Nonlinear cavity enhancement factor showing signif- 
icant increase over the range of m. 
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FIG. 2: Linear inter-disk coupling coefficient for values of m 
ranging from m = 10 (lower line) to m = 50 (upper line) 



enhancement parameter 



r A f = x r 2 \a m /pm^ 



ip(x { ^n d p) 



\ip{x%]n d p)\ 2 pdp. 



( 36 ) 

Both parameters Km an d Tm are in general complex- 
valued, but for high-Q modes their imaginary parts are 
much smaller than 7m, which describes the main contri- 
bution to the radiative losses of our systems. Therefore, 
below we assume that XmiJiM) = 2>?"i(r M ) = 0. The 
resulting equations give a complete ab initio description 
of the nonlinear dynamics of coupled disks in the reso- 
nant approximation. A resonator-induced enhancement 
of nonlinear effects manifests itself through parameter 
r m , which is found to increase by six orders of magnitude 
when the azimuthal number m changes from m = 5 to 
m = 40 (see FigureUJ). This enhancement reflects drastic 
concentration of the field of the WGM modes within its 
effective volume. Fig. [5] shows the linear inter-disk cou- 
pling parameter as a function of the inter-disk distance 
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FIG. 3: Imbalance intensity J± (t) for the two disks sys- 
tem. We have used some generic initial conditions such that 
JS^I 2 = 1.5, |S a) ! 2 = 0.5, |Sf | 2 = |Si 2) | 2 = 0. Plot (a) 
is obtained with x — 2 (self-trapping regime), while plot (b) 
correspond to x = 1 (Josephson oscillation regime). 



for values of the mode order to ranging from to — 10 to 
to = 50. It is interesting, that while the nonlinear pa- 
rameter shows fast growth with to, the linear inter-disk 
coupling actually decreases with to while also becoming 
more short-ranged. This is explained by the fact that 
this parameter is a product of two factors: the Hankel 
function, which grows very fast with to and the radiative 

(r) 

decays rate ^y M which significantly diminishes with to. 
The result presented in this plot shows the interplay of 
these two opposite tendencies. 



III. SELF-TRAPPING TRANSITION IN THE 
DOUBLE DISK SYSTEM 



Non-linear discrete systems, exhibit a transition from 
a Josephson-like oscillatory motion to a self-trapping be- 
havior when the non-linearity strength increases beyond 
some critical value. The phenomenon has been studied 
extensively for a variety of systems ranging from cou- 
pled non- linear waveguides to BEC's in optical lattices 
and biological systems [17H26I ] while recently it has been 
also observed experimentally in the frame of BEC's 24 1. 
Here, we will show that coupled optical micro-discs can 
be used as a prototype system to analyze and observe 
experimentally such phase-transition. 

Our analysis of the temporal behavior of the electric 
field relies on Eq. ([55]) where (to the first approxima- 
tion) we neglect radiative losses (70 = 0). An analytical 
treatment of the dynamics can be achieved if we re- write 
Eqs. (f3"5")) in terms of the Stokes parameters defined as 



Ji = 

j£ s) = 
(y) _ 



'fr I 2 



»(sS(4 2 i)" 



(s£i)*4 2 i) 



'ill 2 



l^l 2 



(37) 



The first of these parameters, Jj., is associated with 
the total field intensity distributed between the coupled 
counter-clockwise (clockwise) and clockwise (counter- 
clockwise) modes of discs 1 and 2 respectively. One 



can show that dJ. 



(0) 



= 0, i.e. 



4 



is a conserved 

0) 



>±'/dr 

quantity. Of special interest is the J±' component 
which describes the intensity imbalance between counter- 
clockwise (clockwise) and clockwise (counter-clockwise) 
modes of discs 1 and 2. Using the Stokes variables Eq. 
(|371) . we can re- write Eqs. ([33)) in the following form: 



cU 



± 



where we 



(It 

introduced 



= J± x B, (38) 
the Stokes vector 3± = 



the pseudo- magnetic field B = 

(2,0,xJ± } + 2xJ^ } ) with X = (V M /h M and redefined 
the time variable as r — > r / Km ■ 

( z) 

For the intensity imbalance between the disks J± , one 
can derive the following equation: 



2 t(*) 



dT 2 



2J^)=Q. 



(39) 



Direct integration of Eq. (|39j) for various initial condi- 
tions indicates that there is always a critical value Xcr 

of 

(z) 

the non-linearity, for which Jj_ shows a transition from a 
regime where it is always positive or negative to a regime 
that oscillates around zero. An example of such behavior 
is shown in Fig. [3£i,b. The former case is associated to 
the self-trapping regime, while the latter one is associated 
with the non-linear Josephson-like oscillations pj], I21I - 
[24j . Obviously, the critical value of the non-linearity Xcr 
above which energy transfer from one disc to another can 
take place, depends on the initial preparation of the field 
excitation in the coupled micro-disc system. 

The smallest value of % cr is realized for initial condi- 
tions, J± \o) = 1 with the norm J±(0) = 1. These corre- 
spond to an initial excitation of the clockwise and coun- 
terclockwise modes of one of the disks. From Eq. (|3"8"|) we 
find that 



J 



(x) L.jW = _|((^*) + j<*))2 + 2j(f)jW) + ^. (40) 



Substituting the first term on the r.h.s. of Eq. (|40l) . 
from Eq. (|3U)) . and using the fact that and are 
symmetric, we eventually get 



d 2 J { + z) 

dT 2 



(d _ 9 x\ T (z) 9 X 2 {z f 
(4 —)J+ + — J+ 



0. 



Equation (HT|) admits the following solution 



(r) 



cn(2r;?7); 77 < 1 
dn(3xr/2;7?- 1 ); 77 > 1. 



(41) 



(42) 



where cn(u, 77) and dn(it, rf) are Jacobian elliptic func- 
tions, and 77 = 3x/4 is the modulus of the elliptic func- 
tion. The value 77 = 1, corresponding to Xcr — 4/3, 



8 



marks a transition from an oscillatory to a self-trapped 
behavior. 

While the previous theoretical analysis allow us to cal- 
culate quantitatively the solutions of Eq. Q38p and derive 
an expression for the critical non-linearity strength, we 
find it useful to provide also a qualitative argument ex- 
plaining the existence of Josephson to self-trapping tran- 
sition. The main observation is based on the fact that 
an initial preparation will be redistributed in a way that 
it will minimize the Hamiltonian function (energy) asso- 
ciated with our system. The latter can be derived from 
the equations of motion (|35[) and has the form 

i 

+ f(|S;| 4 + |S-| 4 ) + 2x(|S;| 2 |Sl| 2 )), (43) 

In the self trapping regime, the energy is concentrated at 
the specifics modes which are initially populated (initial 
field distribution). For such field configuration we get 
%ST = 3%. On the other hand, in the limit of Joseph- 
son oscillations the energy is distributed (on the average) 
uniformly over the whole system. If one assumes that 
S± oc 1/ V2, we get the corresponding value of the en- 
ergy function Hjo — (3x+4)/2. The critical nonlinearity 
Xcr for which the transition from one regime to another 
occurs can be evaluated by equating the two energy dis- 
tributions. Our simple argument gives % cr = 4/3 which 
coincides with the result obtained from the rigorous so- 
lution presented above. 

This analysis, can be further extended to the case when 
the counter-propagating modes in each disk are allowed 
to interact with each other due to, for instance, surface 
roughness induced scattering [U E|. This interaction 
can be described by adding an extra term vmS^ , to the 
left side of Eq. ([35l) . where vm describe the strength of the 
intra-disc coupling at each disc. Using the same analy- 
sis as above, one can estimate the critical value of non- 
linearity. For example, for the initial conditions that we 
have used j'f' = J_ = 1, we find Xcr = 4/3, indepen- 
dently from the value of the intra-disc coupling. Again, 
our numerical calculations agrees with the prediction of 
the heuristic argument. 



IV. LEAKING DYNAMICS 

A natural question is how the self-trapping phe- 
nomenon affects the relaxation dynamics of this system 
once a leakage is introduced in one of the disks. Experi- 
mentally this can be achieved by coupling one of the disks 
to a tapered fiber as discussed in the introduction. This 
will result in the imbalance between decay rates of the 
two disks so that the smaller intrinsic decay rate in the 
second disk can be neglected. Theoretically, we describe 
this situation by restoring a radiative decay rate, 70, only 
in the first of the Eq. p)]) . 




FIG. 4: Temporal behaviour of the total intensity remaining 
inside the two discs when the first one is attached to a fiber. 
Plots (a) and (b) present the rescaled total intensity P(t) 
for initial excitation at the second disk and at the first disk 
respectively. In (a) the leaking constant is 7 = 0.7 while in 
(b) we have used 7 = 0.05. In both cases, we have used 
various non-linearity strengths above and below the critical 
non-linearity parameter. 



The object of interest is the total intensity P(t) = 
jf + </i°\ which is no longer a conserving quantity. 
In order to eliminate effects due to trivial exponential 
decay of the intensity, we present the numerical results 
for the relaxation dynamics in Fig. 01 in terms of the 
rescaled parameter P = Pcxp(7or). Two different relax- 
ation regimes are observed in the short time limit. While 
for x < 4/3 the rescaled norm oscillates with a period 
similar to that of the corresponding closed system in the 
Josephson regime, for x > 4/3 there is an exponential de- 
cay, which is not eliminated by the rescaling procedure. 
The origin of this discrepancy is associated with the self- 
trapping phenomenon. Indeed, if initial excitation is cre- 
ated in the leaky disk and x > 4/3, it becomes trapped 
and starts leaking out at a faster rate. When, however, 
the total intensity decreases and the system moves to 
the Josephson regime the decay process slows down since 
the light intensity oscillates between the leaky and the 
non-leaky disk. On the other hand, if we excite the non- 
leaking disk, then for x > 4/3, the initial decay is slower, 
which manifests itself in the initial growth of the rescaled 
intensity P. 

To achieve a better understanding of the initial decay 
in the self-trapped regime, we have analyzed the eigen- 
value problem associated with Eq. (|35p in the presence 
of the dissipation i.e. 



-ryo<5i,i<I> 



±1 



4> 



0) 



(44) 



(i) |2 



By multiplying each set of the equations in P5j) by the 

Si) * 

corresponding ^> ±1 , and adding the resulting expressions 
for the propagating and counter-propagating modes as- 
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FIG. 5: Decay rates extracted from the analysis of the 
rescaled total intensity shown in Figs. [4^,b. Plot (a) corre- 
spond to initial excitation at the second disk. Plot (b) reports 
the initial decay rates when the excitation is prepared at the 
first disk. In both cases, we have used various non-linearity 
strengths x an d leaking constants 7. The black lines represent 
the theoretical predictions Eq. (|47|) . 



sociated with the two disks, we get: 

£ = +^T*S) -x(i*Si 4 +i*^i 4 ) 

-2 X (|<I>« | 2 |$« | 2 + l^l 2 !^]! 2 ) - Hol^J 2 (45) 

where we have used the normalization |$^ 1 | 2 +|$'f \\ 2 = 1 
for i j and i,j = 1,2. Equating real and imaginary 
parts of both sides of the above equations, we get that 

Xm£ = -7o|$+]| 2 = -7o|$ < L 1 i| 2 . Using the last equality 
together with the normalization condition we can con- 
clude that l^il 2 = l^+il 2 . Therefore, the four sets of 
equations in Eq. (|45j) reduce to the following eigenvalue 
problem 

£#g = -iTb*si-*?l-3 X i*a a *g 

£*<?l = -*f]-3x|$L 2 l| a *L a l (46) 

which can be solve exactly. Specifically, we find that 
for x < X* = (1/3)1/4^ To there are two (equidis- 
tributed among the disks) leaking modes corresponding 
to an energy distributed equally between the two discs 
|$^| 2 = |$^| 2 The associated complex eigenval- 
ues £ have an imaginary part Tm{£) = —70/2 which 
dictates their decay rate. For x > X* two new (non- 
equidistributed) modes appear with 

Xm{£^) = ( 70 /2)(-l ± V /l-(4/3)/[ X 2 + 7o 2 ])- (47) 



These solutions correspond to an non-equal energy oc- 
cupations |$5]| a = |(1 T V 1 " ( 4 /3)/[x 2 + 7 2 ])- As ex- 
pected, the Xm(£( + )) ilm{£^')) decay rate corresponds 
to the mode, which has most of the intensity concentrated 
in the (non-)leaky disk. These rates are in agreement 
with those extracted from our simulations, see Fig. El 
V. CONCLUSION 



Using ab initio approach, we derived the equations that 
describe the dynamics of amplitudes of high-Q WGMs 
in a system of two evanescently coupled microdisk res- 
onators. Taking into account linear coupling only be- 
tween counter-propagating modes in adjacent disks, we 
studied manifestations of self-trapping phenomenon in 
the relaxation dynamics. It should be noted that similar 
effect was also discussed in Ref . [13, Gil within the frame- 
work of quantum electrodynamics. An important distinc- 
tion between the system studied in this paper and those 
considered in Ref. [Hd 

lies in the mechanism responsi- 
ble for enhancement of nonlinear effects. In Ref. [T3, Gil 
this enhancement was due to resonant interaction of pho- 
tons with material excitations such as atoms or excitons, 
while in this work the enhancement occurs due to small 
volume of WGMs. As a result, the nonlinear effects con- 
sidered in this paper can be observed at room tempera- 
ture. For instance, we found that in the presence of the 
decay rate imbalance between the coupled disks, there 
exists an anomalous relaxation behavior similar to the 
one discussed in Ref. [27l - l3ll |. This result identifies the 
system of coupled optical resonators as a convenient plat- 
form for experimental study of this phenomenon, which 
so far, has eluded experimental observation. 
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